

//  Include "petscksp.h" so that we can use KSP solvers.  Note that this file
//  automatically includes:
//     petsc.h       - base PETSc routines   petscvec.h - vectors
//     petscsys.h    - system routines       petscmat.h - matrices
//     petscis.h     - index sets            petscksp.h - Krylov subspace methods
//     petscviewer.h - viewers               petscpc.h  - preconditioners
//
#include "petscksp.h"


#include ".\parallelFE.h"

//
/////////////////////////////////////////////////////////////////////
parallelFE::parallelFE(void)
{
  int ierr;

  PetscInitialize(NULL,NULL,(char *)0,NULL);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
  CHKERRQ(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);
  CHKERRQ(ierr);
}

//
/////////////////////////////////////////////////////////////////////
parallelFE::~parallelFE(void)
{
  int ierr = PetscFinalize();
  CHKERRQ(ierr);
}


// Set dimension of problem
/////////////////////////////////////////////////////////////////////
int parallelFE::beginAssemble (int N)
{
  int ierr;

  u     = (Vec *) malloc (sizeof (Vec));
  b     = (Vec *) malloc (sizeof (Vec));
  ustar = (Vec *) malloc (sizeof (Vec));
  A     = (Mat *) malloc (sizeof (Mat));

  //   Create stiffness matrix
  ierr = MatCreate (PETSC_COMM_WORLD, A);
  CHKERRQ(ierr);
  ierr = MatSetSizes (*A, PETSC_DECIDE, PETSC_DECIDE, N,N);
  CHKERRQ(ierr);
  ierr = MatSetFromOptions(*A);
  CHKERRQ(ierr);

  //   Create right-hand-side and solution vectors
  ierr = VecCreate (PETSC_COMM_WORLD, u);
  CHKERRQ (ierr); 
  ierr = VecSetSizes(*u,PETSC_DECIDE,N);
  CHKERRQ (ierr); 
  ierr = VecSetFromOptions(*u);
  CHKERRQ (ierr);
  ierr = PetscObjectSetName((PetscObject)u, "Approx. Solution");
  CHKERRQ (ierr);
  
  ierr = VecDuplicate (*u, b);
  CHKERRQ (ierr);
  ierr = PetscObjectSetName ((PetscObject)b," Right hand side");
  CHKERRQ (ierr);
  ierr = VecDuplicate (b, &ustar);
  CHKERRQ (ierr);
  ierr = VecSet (u, 0.0);
  CHKERRQ(ierr);
  ierr = VecSet (b, 0.0);
  CHKERRQ (ierr);

  return 1;
}


// insertLocalStiffnessMatrix
/////////////////////////////////////////////////////////////////////
int parallelFE::insertLocalStiffnessMatrix (int ndof, int *idx, double *Ke)
{
  int ierr = MatSetValues (*A, ndof, idx, ndof, idx, Ke, ADD_VALUES);
  
  CHKERRQ(ierr);

  return 1;
}


// endAssemble
/////////////////////////////////////////////////////////////////////
int parallelFE::endAssemble (void)
{
  int ierr;

  ierr = MatAssemblyBegin (*A,MAT_FINAL_ASSEMBLY);
  CHKERRQ (ierr);
  ierr = MatAssemblyEnd (*A,MAT_FINAL_ASSEMBLY);
  CHKERRQ (ierr);

  ierr = VecAssemblyBegin (b);
  CHKERRQ (ierr);
  ierr = VecAssemblyEnd (b);
  CHKERRQ (ierr);
}
